Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.
In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv
Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates. Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.
3) Add sample random effect structured by cohort to allow cohort effects on growth.
18.1 By Cohort, modelNum 3: phi(i,t) * g(i,t), p(i,t) model
Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)
Model code is in ./models/production/modelProduction_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/modelProduction_OB_makeFile.R
Functions for this qmd file in ./models/qmdProduction_OB_functionsToSource.R
Model runs:
18.1.1 How many ageInSamples to include?
Code
all <-tar_read(cdWB_electro_target)table(all %>%filter(river =="wb obear")|> dplyr::select(ageInSamples))
library(targets)# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html# # # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAssource('./models/production/modelProduction_OB_functionsToSource.R')source('./models/production/qmdProduction_OB_functionsToSource.R')modelNum <-3# phi * growth
18.1.3 Model code
Code
# all cohorts have the same code. show hereload('./models/production/runsOut/production_OB_model_3_2002__mostRecent.RData')d$modelCode
{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]],
sd = 2)
weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = weightSD)
for (t in (first[i]):(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]/90
weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1],
sd = weightSD)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ dnorm(0, sd = 1000)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i,
t] * weight[i, t] + phiBeta[2, i, t] * weight[i,
t] * weight[i, t]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
for (b in 1:2) {
phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
}
}
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
phiIntTOut[t] <- ilogit(phiIntT[t])
phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
weightZ01[i, t] <- weight[i, t] * z[i, t]
}
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}
#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same ehojs_define(summary_OB_OJS_all = sAll)#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))
Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
18.1.8 Plot survival
Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
# Production model - river 4 (IS)Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.In all the models below, 1 = not observed and 2 = observed in the input *encounter data*.\Encounter data are available [here](https://drive.google.com/drive/folders/1UmPv49xL-mzUBndqjRISylNO3q3Zs-cm) in the `eh.csv` file. Weight data are in `weight.csv````{r}#| label: librariesModelsNimbleProduction#| echo: falselibrary(getWBData)library(lubridate)library(kableExtra)library(GGally)library(nimble)library(nimbleEcology)library(MCMCvis)library(tidyverse)library(targets)library(shiny)library(tibble)``````{r}#| label: globalModelsNimbleProduction#| include: falseknitr::opts_chunk$set(warning =FALSE, message =FALSE)```Using model #3 from `modelsCMR_Growth_NN_OB.qmd` as a staring point for the models, but adapting the model by\Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates. Prob not 2) Loop over first\[i\]:lastAIS so fish have the chance to survive to large size.\3) Add sample random effect structured by cohort to allow cohort effects on growth.## By Cohort, modelNum 3: phi(i,t) \* g(i,t), p(i,t) modelModel with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)Probability of survival (phi) model structure:``` logit(phi[i,t]) <- phiInt[i,t] + phiBeta[1,i,t] * weight[i,t] + phiBeta[2,i,t] * weight[i,t] * weight[i,t] + phiBetaCohort[cohort[i]] ```Probability of capture (p) model structure:p(t,i) \<- pInt(t,i)Growth rate (gr) model structure:gr(t,i) \<- grInt(t) + cohort(i)Survival/growth rate interaction:MultiplicativeModel code is in `./models/production/modelProduction_OB_functionsToSource.R`\Model is run 'by hand' in `./models/modelProduction_OB_makeFile.R`\Functions for this qmd file in `./models/qmdProduction_OB_functionsToSource.R`\Model runs:\### How many ageInSamples to include?```{r}#| label: productionall <-tar_read(cdWB_electro_target)table(all %>%filter(river =="wb obear")|> dplyr::select(ageInSamples))cohorts <-2002:2014```### Retrieve model results```{r}#| label: productionSource#| cache: falselibrary(targets)# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html# # # To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAssource('./models/production/modelProduction_OB_functionsToSource.R')source('./models/production/qmdProduction_OB_functionsToSource.R')modelNum <-3# phi * growth```### Model code```{r}# all cohorts have the same code. show hereload('./models/production/runsOut/production_OB_model_3_2002__mostRecent.RData')d$modelCode``````{r}#| label: productionFunctiongetSummaryByCohort <-function(cohort, modelNum =3) {# Load the dataload(paste0('./models/production/runsOut/production_OB_model', '_', modelNum, '_', cohort, '__mostRecent.RData')) # Input data eh <-tar_read_raw(paste0('eh_OB_', cohort, '_target'))get_MCMCTrace_OB(d)MCMCplot(object = d$mcmc, params =c("phiIntTOut"), ref_ovl =TRUE)MCMCplot(object = d$mcmc, params =c("phiBetaT"), ref_ovl =TRUE)MCMCplot(object = d$mcmc, params =c("p"), ref_ovl =TRUE)MCMCplot(object = d$mcmc, params =c("grIntT"), ref_ovl =TRUE)# Create the summary list summary_list <-list(eh = eh,stats =kable(as_tibble(d$runData), caption ="Run statistics"),runTime =paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")),#traces = get_MCMCTrace_OB(d),summaryMod3_tt_growth =MCMCsummary(object = d$mcmc, params =c("phiIntTOut", "p", "grIntT", "phiBetaT"), round =3) %>%rownames_to_column(var ="var"),summary_OB =get_summary_OB(d, eh) )return(summary_list)}#s=summary_OB[, order(summary_OB$ind, summary_OB$t)]```### Summaries by cohort::: {.panel-tabset}#### 2002```{r}#| label: s2002#| cache: TRUEs2002 =getSummaryByCohort(2002, modelNum)```#### 2003```{r}#| label: s2003#| cache: TRUEs2003 =getSummaryByCohort(2003, modelNum)```#### 2004```{r}#| label: s2004#| cache: TRUEs2004 =getSummaryByCohort(2004, modelNum)```#### 2005```{r}#| label: s2005#| cache: TRUEs2005 =getSummaryByCohort(2005, modelNum)```#### 2006```{r}#| label: s2006#| cache: TRUEs2006 =getSummaryByCohort(2006, modelNum)```#### 2007```{r}#| label: s2007#| cache: TRUEs2007 =getSummaryByCohort(2007, modelNum)```#### 2008```{r}#| label: s2008#| cache: TRUEs2008 =getSummaryByCohort(2008, modelNum)```#### 2009```{r}#| label: s2009#| cache: TRUEs2009 =getSummaryByCohort(2009, modelNum)```#### 2010```{r}#| label: s2010#| cache: TRUEs2010 =getSummaryByCohort(2010, modelNum)```#### 2011```{r}#| label: s2011#| cache: TRUEs2011 =getSummaryByCohort(2011, modelNum)```#### 2012```{r}#| label: s2012#| cache: TRUEs2012 =getSummaryByCohort(2012, modelNum)```#### 2013```{r}#| label: s2013#| cache: TRUEs2013 =getSummaryByCohort(2013, modelNum)```#### 2014```{r}#| label: s2014#| cache: TRUEs2014 =getSummaryByCohort(2014, modelNum)```:::### Combine cohort summaries```{r}#| label: bindSummariessAll <-tibble(.rows =0) |>add_column(!!!s2002$summary_OB[0,])# Fill tibblefor (cohort in cohorts) { summary_obj <-get(paste0("s", cohort))$summary_OB summary_obj$cohort <- cohort sAll <-bind_rows(sAll, summary_obj)}sAll$cohort_ind <-paste0(sAll$cohort, "_", sAll$ind)``````{r}#| label: ojsDefine#| cache: FALSE#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same ehojs_define(summary_OB_OJS_all = sAll)#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))```### Select cohort(s):```{ojs}//| label: summary_OB_OJS_all_Tsummary_OB_OJS_all_T = transpose(summary_OB_OJS_all)``````{ojs}//| label: selectCohort//| import { range } from "d3";viewof selected_cohort = Inputs.select(d3.range(2002, 2015), { label: "Which cohort?", value: 2002, step: 1, multiple: true})``````{ojs}//| label: filterCohortsummary_OB_OJS_all_T_selected_cohort = summary_OB_OJS_all_T.filter((d) => (selected_cohort.includes(d.cohort)))```### Select one or more individuals:```{ojs}//| label: selectInduniqueInds = [...new Set(summary_OB_OJS_all_T_selected_cohort.map(d => d.cohort_ind))].sort();viewof selectInd_mod3 = Inputs.select(uniqueInds, { label: "Which fish?", value: 1, step: 1, multiple: true})summary_OB_OJS_all_T_selected_cohort_ind = summary_OB_OJS_all_T_selected_cohort.filter((d) => (selectInd_mod3.includes(d.cohort_ind)))```Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort.### Plot survivalBlack dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish's cohort.```{ojs}//| label: plotZ Plot.plot({ width: width, height: 350, inset: 10, color: { scheme: "greys" }, x: { label: "Age/season combination" }, y: { label: "Probability of survival" }, marks: [ Plot.frame(), Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "pSurv" }), Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "enc8", fill: "orange" }), Plot.line(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "pSurv" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "first", y: 1, "stroke": "green" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "last", y: 1, "stroke": "red" }) , Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, Plot.selectFirst({ x: 11, y: 1, text: d => d.cohort }) ) ], facet: { data: summary_OB_OJS_all_T_selected_cohort_ind, x: "ind" } })``````{ojs}//| label: plotWeightPlot.plot({ width: width, height: 350, inset: 10, color: { scheme: "greys" }, x: { label: "Age/season combination" }, y: { label: "Standardized body mass (mg)" }, marks: [ Plot.frame(), Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "mean_weight" }), Plot.line(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "mean_weight" }), <!-- Plot.dot(d, { --> <!-- x: "t", --> <!-- y: "mean_gr", --> <!-- fill: "green" --> <!-- }), --> Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind, { x: "t", y: "weightDATA_std", fill: "orange" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "first", y: 3, "stroke": "green" }), Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind, { x: "last", y: 3, "stroke": "red" }), Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, Plot.selectFirst({ x: 1, y: 4, frameAnchor: "top-left", text: (d) => "Cohort = " + d.cohort }) ), Plot.text(summary_OB_OJS_all_T_selected_cohort_ind, Plot.selectFirst({ x: 1, y: 3.5, frameAnchor: "top-left", text: (d) => "Residual = " + d.meanResid }) ) ], facet: { data: summary_OB_OJS_all_T_selected_cohort_ind, x: "ind" } })```